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AN ITERATIVE METHOD FOR FINDING THE 
TRANSIENT RESPONSE OF AN AUTOMATICALLY 
COI^TROLLED FISSILE HAVING NON-LINEAR STIFFNESS 



SUMMARY 



An equation of missile motion in the pitching plane was found, 
assuming a non-linearity in aerodynamic stiffness in the form of a 
cubic term. An iterative method for approximating the angle-of- 
attack response to a step displacement of elevator angle was 
formulated. The results of applications of this method were 
compared with computer solutions; convergence of the approximations 
wopI demonstrated. 

A simple feedback control system was designed, and the iterative 
method was extended to the case of a control system step input. 

The method was thought to have limited practical value because 
of the labor required in using it. 
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AN ITERATIVE METHOD FOR FINDING THE 



TRANSIENT RESPONSE OF AM AUTOMAT ICALLY 
CONTROLLED MISSILE HAVING NON-LINEAR STIFFNESS 



INTRODUCTION 



Important non-linearities which are often present in low aspect 
ratio configurations are those of the lift versus incidence and the 
pitching moment versus incidence characteristics. In missile design 
these non-linearities, as well as any others present, are likely to 
be ignored initially, simply because linear system design is well 
developed and relatively easy. Automatic computing techniques may 
then be used to determine the behaviour of the system when the non- 
linearities are included. 

It would be a great help for the designer to be able to find 
the transient response of a contemplated design without going to the 
computer. It was for this reason that the present investigation 
was undertaken. By expressing the pitching moment as the sum of a 
linear tern and a cubic term in angle-of-attack, ignoring the 
non-linearity in lift because it is likely to be less important, 
an equation amenable to iterative solution can be formulated. The 
first approximation to the transient response may be expressed as a 
Fourier series, which must be cubed to perform the first iteration. 
The main effort of the development will be to formulate a practical 
method for performing the cubing process, for the first and 
subsequent iterations, and to systematise the steps necessary for 
carrying out the iterations. 

First, an equation relating angle-of-attack to elevator angle 
will be derived, and the response to a step displacement of elevator 
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angle will be found. Later, a simple automatic feedback control 
system will be designed for the purpose of formulating a method 
of finding the response to a control system input. 

The author takes this opportunity to express his gratitude for 
the generous assistance given by Mr. P. A. T. Christopher, who 
instigated and supervised this investigation. 
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EMGLISH-LbTiER N0I1EIJCLA.TURE 

coefficient of restoring term in Sec. 3 
coefficient of damping term in Sec. 3 
coefficients of terms in the cube of Fourier series 
in Sec. 2 

gain of frequency response at frequency n>^ F 

coefficient of cubic term in Sec. 3 

coefficients in the Fourier series for 

phase shift of frequency response at frequency n*) F 

coefficients in Sec. 2 

phase angles in the Fourier series for \ ^ 
l/K^ times the transfer function of the compensator- 
amplifier 
lift coefficient 
pitching moment coefficient 

partial derivative of My with respect to ^ ^ 

coefficient of vane damping term 

decibels 

d/dt 

missile body diameter 
actuating error 

coefficients in trigonometric identities 
phase angles in trigonometric identities 
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H(s) 



i 



transfer function used in closed-loop iteration 

f(-i) 



K i 






moment of inertia of missile about y axis 
moment of inertia of vane about its c*g. 
zero frequency gain of the aerodynamic transfer 
function 

zero frequency gain of compensator-amplifier 
gain of input amplifier 

distance between missile c*g. and vane c*g* 
lift 



M 

fc'(w) 



M 



C<* 



M, 

n 

P(s) 

q 

Q 



mass of missile 

aerodynamic pitching moment on missile about y axis 
aerodynamic pitching moment on vane about its c.g. 
component of moment due to w 
component of moment due to 

coefficient of the linear term in the cubic expression 
for pitching moment 

coefficient of the cubic term in the cubic expression 
for pitching moment 

partial derivative of M with respect to ^ 
harmonic 

l/K^ times aerodynamic transfer function 
angular velocity of missile about the y axis 

step function in Sec* 3 
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Q stagnation pressure in appendices I and VI 
IQ ( magnitude of step function 
non-dimensional step input 
|Q | magnitude of non-dimensional step input 

R ratio of non-linear to linear steady-state values 
s Laplace transform variable 
S wing area 
t time 

t^ time measured from start of second half cycle 
T period of square wave 

T(d) differential operator in equation for missile 
aerodynamics 

u component of missile c.g. velocity in x direction 
U(D) differential operator in equation for missile 
aerodynamics 
V velocity of missile 
V(s) vane transfer function 

w component of missile c.g. velocity in z direction 
W(s) actuator transfer function 

Z aerodynamic force component in the z direction 

Z partial derivative of Z with respect to w 

Z^ partial derivative of Z with respect to 

Z^ partial derivative of Z with respect to 
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GREEK- LETTER NOMENCLATURE 



missile angle-of-attack 

demanded missile angle-of-attack 

magnitude of input step function 

reference input 

vane angle-of-attack 

approximations to 

corrections to approximations 

response of linear system to unit step input 

coefficient of the cubic term 

angle between missile axis and vane axis 

damping ratio 

damping ratio of missile 

damping ratio of vane 

damping ratio of elevator actuator 

elevator deflection angle 

magnitude of step elevator input 

demanded elevator angle 

angle between missile axis and space reference 
angle between vane axis and space reference 
time constant in the missile aerodynamics transfer 
function 

fundamental frequency of square wave 
natural frequency of linear system 
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uj>^ natural frequency of linear missile 

a) natural frequency of vane 

** L 

natural frequency of elevator actuator 



SUBSCRIPTS NOT SHOW ELSEWHERE 



a evaluated function 

b function to be derived from an evaluated function 
L linear 
N non-linear 
S steady-state 



8 



SECTION 1: THE EQUATION OF MISSILE MOTION IN PITCH 

The system of axes and sign conventions to be used are 
illustrated in Pig. 1. Body axes are used. The missile to be 
considered lias a fixed wing and an all-moving tail. It may or may 
not be of cruciform conf iguration, since pure pitching motion will 
be assumed. 

The rigid body equations of motion in the z direction and about 
the y axis (which passes through the center of mass) are, 
respectively, 



Z = m(w - qu) 


(1-1) 


M = Iq 


(1-2) 



The same equations would, of course, arise if the general three 
dimensional equations of motion were written down and the confining 
assumptions of the present motion imposed. 

The component in the z direction of the total force on the 
missile is taken as 

z = V + z i n d-3) 

The simplifying assumption has been made here that the gravitational 
force may be neglected because it is small compared with the 
aerodynamic forces available. The contribution to Z due to q is 
omitted on the assumption that it is small. Its inclusion would 
not change the nature of the analysis, but it would be dropped 
eventually anyway because values for it, for the conf iguration and 
flight conditions considered, are not available. Also implicit in 
Eq. 1-3 is the assumption that the malalignment force is zero, 
i.e. Z is zero when w and 



are zero. 
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The aerodynamic moment is taken as 

M = M'(v)+M n (1-4) 

where M'(w) is to be read, "The component of moment due to w." 

This is the non-linear term. The same remarks as those following 
Eq. 1-3 apply to the exclusion of components of moment due to q 
and to malalignment. 

Now we take 



u = V (1-5) 

which is valid for small angles of attack, and we assume that 
V is constant. Then, for small angles of attack, 

< -7 d-6) 

and 

z„ - l* d-7) 

W V 

Substituting Eqs. 1-5 and 1-6 in Eq. 1-1 gives 

Z = m(V<* - qV) (1-8) 

Substituting Eqs. 1-6 and 1-7 in Eq. 1-5 gives 

z = z^*+z ^ (1-9) 

Also we replace M'(w) by M"(«<), where M"(o0 is to be read, 

"The component of moment due to &C . 11 Thus Eq. 1-4 becomes 

M = M"(®<) + n (1-10) 

Combining Eqs. 1-8 and 1-9 gives 

Z* ^ + Z ^ H = m(V<k - qV) (1-11 ) 

Combining Eqs. 1-2 and 1-10 gives 

M"(o<) + ^ = Iq 



( 1 - 12 ) 
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Eliminating q between Eqs. 1-11 and 1-12 gives 

mlV e< - IZ^ ^ - MVM"(®<) = IZ >\ + mVM^ ^ (l-13) 

Taking 

M"( «*) = <* + A* (1-14) 

equation 1-13 becomes 

mlV ®< - IZ^ - mVM^ <K. - mVM^i = IZ^ ^ + mVM^ ^ (l-15) 

Dividing Eq. 1-15 through by mlV gives 



©< - <X - 

mV 


* * % 1 . 

o< “ = h + 

I I mV 1 


1 

- 1 


(1-16) 


The constants , K,, S' 


and /3 are 


now defined 


as follows: 








II 


_ 

i 




(1-17) 


P = 


. Z * 

2mVwl^ ( 




(1-18) 


II 

i — 1 
& 






(1-19) 


I 10 Z 




z = 


Z h 




(1-20) 




mVK^ 




ft = 


_ 




(1-21) 



I 



Making these substitutions in Eq. 1-16 gives 

o < f 2. + L O"' * +/3 <*“* - K' i + Z ( 1 - 22 ) 

In Ref, 1, E. G. Brown-Edwards calculated the aerodynamic 
coefficients for a typical missile configuration. This missile 
is sketched in Pig, 2. In Appendix I Brown Edwards 1 results are 
given and are used as the basis of the calculation of the constants 
in Eq, 1-22. The body diameter was taken as one foot and the flight 
Mach number as three. The calculations for a flight altitude of 
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50,000' are shown in the appendix. The results are tabulated in 
Table I. This table also contains the results for sea level 
flight, the calculations for which are ncrt shown. 

SECTION 2 : THE TRANSIENT RESPONSE OF THE MISSILE TO A 

STEP DISPLACEMENT OF ELEVATOR ANGLE 



Equation 1-22 relates the angle-of-attack °< to the elevator 
angle ^ . We wish to find the variation of °< with time upon 
the application of a step displacement in ^ , initial values all 
being zero. No analytical solution for this non-linear problem 
exists; an iterative method will be sought. 

Equation 1-22 is, in operational fora, 



(D 2 + 2 5^ D + < )o< +/3o< > = K]l uJ^, 2 (l + ^ D) 

This may be written: 

T(D) C< + ft •<* = u(D)*i (2-2) 

where 

T(D) = D 2 + 2 ^ D + uC, ( (2-3) 

U(D) = * 1 ' A C t (1+ H) (2-4) 

Now an approximation to the solution of Eq. 2-2 is the 
solution (t) of 

T(D) •<, = U(D) ^ (2-5) 



This equation is linear, and could be solved analytically, however 

we will use an approximate method. In Ref. 1, where an attempt is 
made to solve the same non-linear problem, Eq. 2-5 .is solved 
analytically. When this is done the ensuing steps become quite 
tedious. We hope to find an advantage in taking another course. 
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We find o< ( (t) in the form of a Fourier series using the 
approximate method given by Wass and Haymen in Ref. 2 and described 
in Appendix II of this thesis. The method requires that the 

frequency response of Eq. 2-5 be determined. The reader who is 
unfamiliar with frequency response technique will find it explained 
in Section 5> in connection with the design of the control system. 
The Wass, Hayman method gives (t), 
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approximately, in the form, 

<(*> - 111 [*%* k ™ ^ t - V] ( 2 - 6 ) 

/n. = 1> 3, ... ,11. 

where <*) , is a known constant, hi is the magnitude of the step 
input, and the A’s and B's are known constants - gains and phase 
shifts of the frequency response. If positive damping exists, 
this response is a damped oscillation, with a steady-state value of 

(oMs-InK ( 2 -7) 

The true steady-state value ( ) may be calculated from Eq. 1-22 

by setting all the terms involving derivatives equal to zero: 




(*) S 




a. 

= K. 

1 



hi 



(2-8) 



The solution for ( °< )^ is easy, though not analytic. 
Now if we take, as a new approximation to o<(t), 



<* 2 (t) 



o< (t) 

»niA 0 



= ilk \ * 2 T... 

a I 2 

o L 

then its steady-state value is 



i 



(2-9) 



o< ) 

2 J S 




(“Os 



o 

Since °C,(t) has the true steady-state value, there is some 
reason for thinking it is an improved approximation for c<(t). 

Some further justification is presented in Fig, 3- The curves in 
this figure represent the solutions, by digital computer, of an 
equation essentially the same as (l-22). The equation is discussed 
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in detail in Section 3* The solutions for the linear case and for 
several non-zero values of the coefficient in the cubic term are 
plotted. In relation to the present example, the linear solution 
corresponds to ^^(t), and the non-linear solutions to oi (t). 

In Fig. 4, each of the non-linear solutions has been scaled up to 
give steady-state values equal to (°<^) - an inversion of the 

process represented by Eq. 2-9- It may be seen that this scaling 
process apparently brings the non-linear solutions into closer 
agreement with the linear one in a general sense, not just the steady- 
state values. Thus ^^(t), a scaled-down version of °^t), will 
be taken as an improved approximation to c<(t). For brevity let 

R = 






so that 



<* 2 (t) = R* x (t) 
For later convenience let 






( t ) = Y A< n sin ( 'n lo t + b ) 

y 2 TrZi 7^. 



(, 2 - 10 ) 



( 2 - 11 ) 



( 2 - 12 ) 



i.e. the response of the linear system to a unit step input, 
ihus, 

°Vt) = hi (t) 



(2-13) 



and 



o< At) = R lul (t) 

* 

Row if /3 is sufficiently small, a closer approximation 
to o( (t) is the solution o<^(t) of 



(2-14) 



T(B>) o< 3 = U(D) rj - (S 



U(D)»j - ft R 3 l^| 3 °< 3 



(2-15 
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Or, if ^^(t) is the solution of 

T (°)°V = -£<*{ = -/* r3 hi 3 

then by the principle of linear superposition, 



( 2 - 16 ) 



°< 3 (t) = ^(t) + ^ 5 *(t) (2-17) 

To solve (2-16), a method will be developed for expressing 
as a Fourier series. 



3 _ 

1 ., 



k f + B n >| 



/Vt =1, 3> . • • , 11. 



The expanded cubic consists of the cube of each of the original 
terms, plus three times the product of the square of each tern and 
each of the other terms, plus six times the product of all 
combinations of three unlike terms. Showing a few terms of the 
expansion and introducing a notation for their coefficients we 
have : 



= a 000 A 0 + a 001 A 0 A 1 sin ^ U> F t + B l^ 



a lll A f 



sm 






V 



+ a 115 k 2 sin 2 (u) p t + B 1 ) sin (3^pt + B^) + 



+ a. A. A .A, sin (iul n t + B.) sin ( j u), ,t + B.) sin (ku) t + B, ) 

ljk l j k ' F l F F k' 



Each term is identified by a subscript i, j, k. The coefficients 
a ijk ^’ a ^ e accoun t of the kind of term (e.g. the multiplier 6 for a 
term which is the product of three unlike terms) and the factors 



( 2 - 18 ) 
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l/2, 2/tt , and l/n. The coefficients a^ k are tabulated in Table II. 
Any product of sine terms can be expressed in the form. 



E o + sin ('Opt + P 1 ) + E 2 sin ( 2 uj^t + F 2 ) 

+ ... + E„ sin (a «J t + F„ ) + . . . 

Por example, in the 1 , 3,5 term we have the identity: 
sin (‘Opt + Bj) sin (3ui p t + B^) sin (5«O p t + Bj.) 

= - i sin (u) p t ) + i sin (3^ p t + Bj-B^B,-) 

+ 7 sin ( 7 <-Opt -B^+B^+B,-) - i; sin (9 uipt +B 2 +B^+B^) 



In this case 



E ]_ = - T r E3 = i , E^ 



JL 

4 f 



= 



E o = E 2 = E 4 - % - E 6 - E 8 - E 10 = 



... = 0 



P 1 “ “ B 1 B 3 + B 5 

P 3 = 

F ? = - B 1 + B ? + B 5 

F g = B ! + B 3 + B 5 



Using this notation a large number of these identities may be 
tabulated concisely, and this has been done, for the combinations 
necessary here, in Table III. 



We adopt the notation#: 
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sin (iuyt + B i ) sin (j u> p t + B^) sin (ku) p t + B^) 

= E / . \ + E, / . .. \ sin (u3 t + F n / . .. \) + . . . 

o(ijk) l(ijk) F l(ijk)' 

Also, we will write, for later convenience, 

sin (iio F t + B i ) sin ( j iO p t + B^) 

= E / . . \ + E, / . . \ sin ( to t + F, / . . \ ) + . . . 

o(oij) 1 (oij ) F l(oijr 

and 

sin (iu> p t + B i ) 

= E / . \ + E .. i . \ sin (to^t + F, / . \) + . . . 

0 ( 001 ) l(ooi) v F l(ooi) / 



(2-19) 



( 2 - 20 ) 



( 2 - 21 ) 



The latter expression is a trivial use of such a notation since the 
right hand side reduces to sin (i*o p t + B^), but its use helps 
provide a systematic tabulation. 

Using the identities represented by Eqs. 2-19 , 2-20, and 2-21 

'Z 

in Eq. 2-18 gives as a series of constants and simple sine 

i o 

terns. Since sine terms of like frequency can be combined into a 
single sine tenn of that frequency, we have: 

* i u 3 = b Q + bi sin (u^t + C x ) + b 2 sin (2<*i p t + C 2 ) (2-22) 

+ . . . + b^ sin ( /vv tO p t + ) + . . . 



where 
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+ a ijk E o(ijk) A i A j \ + •" 

b i sin (“V + C i^ = a ooi E i(ooi) A o A i sin ^F* + F l(00l)^ 
+ a 003 E l(003) A o A 3 Sin ^ U> F t + F l(003)^ + 

+ a ijk E l(iJk) A i A j A k 8in(w, P t+ F l(ijk) } + '** 



b sin (n <jO,,t + C ) = ... 

n F n 

+ a.., E \ A. A. A. sin (n <^> c t + F /. ., \) + ... 

ljk n(ijk) i j He F n(ijk)' 

It is possible to shorten the notation conveniently by 
letting 



^n(ijk) a ijk E n(ijk) (2—23) 

The values of c n (ij] 1 ) rosy be determined from the tabulations of 
a ijk and E n(i,jk)* Accordingly, values of c^^ (or simply c) 

are tabulated in Table IV. The coefficients are all listed as 
positive. The coefficients for which the values of E were negative 
were changed in sign by adding 180° to the corresponding angles F, 
which are also tabulated. It will be seen that like harmonics 
have been grouped together. 

Table IV does not include coefficients for all possible terms. 

The elimination of less important terms was based on the effect of 
each term on °^(t) (the approximation towards which we are working) 
in the numerical example which follows. It had been hoped, at the 
outset of this investigation, that the number of terms needed would 
be small. But it was found necessary to retain a large number of 
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terms, up to and including 9th harmonic terns. The relative 
importance of the smallest terms retained can best be seen in the 
example of this section. 

The fact that the requirements of a particular numerical 
example must be injected at this point is a deficiency in the method, 
since the choice of retained terms is not necessarily the best choice 
for other examples to which this method will be applied. However, 
to follow a more general course would be impractical. 

To obtain values for b^ and for a particular problem, 

the products A^ and the angles are 

calculated using the values of gain A and phase shift B from the 
frequency response curves, then like harmonics are combined. The 
combination of several sine terms in the same frequency may be done 
most easily by graphical vector addition. This is described in 
Appendix III. 

rz 

The Fourier series for ‘K, may be computed with minimum 

1 v 

effort by using the form of Table IV. All the necessary results 
of the preceding development are summarised there, and the required 
computations are indicated by its layout. To use tne table, 
proceed as follows: 

1. Enter the values of gain A and phase shift B of the 
appropriate frequency response. Once this is done, no 
further reference need be made to any information other 
than that contained in Table IV. 



2. Calculate the products A^ A^ A^. 



3. Calculate the products c , . \ A. a. A, . 

^ nujkj i j a 

table cAAA is written for c , . \ A. A. A, ) 

n(ijk) i j k' 



(in the 
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4 . Calculate the angles 

5. Combine like harmonics (e.g. by the method of 
Appendix III), and enter the results, b n and C n , 
in the space provided. 

'Z 

Step 5 above gives the Fourier series for , . Multiplying 

each of the coefficients of the series by the constant - IT hr 

3 3 3 

gives the Fourier series for - Rr hr By Eq. 2-16, 

ot ,*(t) may be found by applying the frequency response ^/T(A.ui). 

The next approximation to ®<(t), if it is desired to go 
further, is °<\ (t), the solution of 

T(B)<*, = U(D) 1 - ft <* 3 3 (2-24) 

*z 

The determination of a Fourier series for may be obviated 

'Z ** 

as follows. Adding and subtracting on the right hand 

side of (2-24) : 

T(B)^ 4 = U(D) ^ -/a^ 2 3 + (3 ^2 - °</) (2-25) 

Let °^*(t) he the solution of 

T(D) = (i (°< 2 3 - ^ 3 3 ) (2-26) 

Then by superposition of the solutions of Eqs. 2-15 and 2-26, 



^(t) = °< 5 (t) + °^ 4 *(t) 



(2-27) 



•z "Z 

The Fourier series for ( °^ ' - ^ . ) need not be determined with 

'Z 

as great an accuracy as would be required for that of ; • To 
solve Eq. 2-26, and ° ^ are each evaluated and cubed, point 
by point. The difference is then plotted, and a numerical harmonic 
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analysis is performed. The frequency response is applied 

to the resulting courier series, giving Instead of 

evaluating ^ 9 (t) and cubing each point, the series which has been -f j 

obtained for °v> could oe evaluated « A be cause cubing a few numbers 

is much easier than evaluating a Fourier series. but the latter 

method has a certain superiority in that if the series obtained for 

o( ^ is not a good representation of the true cube of due 

either to a mistake or to a deficiency in the cubing method, then 

the error will be compensated only if the latter method is used. 

A deficiency in the cubing method would probably exist if the linear 
system frequency response was very different from that used in 
formulating the cubing method. 

Improved approximations, ^ etc., would be found in the 

same way as is 

A method is shown in Section 6, in connection with the control 
system response, for simplifying the calculations when the responses 
to more than one input magnitude are required. An analogous 
method can and should be applied here if more than one input 
magnitude is used. 

We now apply this method for the case of missile flight at sea 
level and an input magnitude of .25 radian. This case was chosen 
because its analog computer solution was available. The constants 
in Eqs.^2-3? and 2-4 are those tabulated for sea level in Table I, 
with one exception: 7 is taken as zero. The effect of that 

term on the transient response was found to be negligible, having 
an influence on the frequency response only at very high frequencies, 
and it was omitted from the analog computer setup. It was 
retained, in the development just shown, for generality. 
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The frequency response of Eq. 2-5, 





T 



(i*>) 2 + 2% J S 1 W (ito) 




1 




(2-28) 



1 - 



is plotted in Fig. 5* Gains are plotted in decibels, and phase 
shifts in degrees. Following the Wass, Hayman method for determining 
linear transient response, we take U> ( to be one fifth of the 
frequency at peak gain. The peak occurs at 



The gains and phase shifts of Eq. 2-28, at frequencies 

nu)p where n = 0, 1, 2, ..., 13, are tabulated in Columns 2 and 
3 of Table V, with values of n in Column 1. The values of gain 
are absolute, not decibels. This procedure is necessary to make 
the computations systematic (this differs slightly from the Wass, 
Hayman procedure). 



oO = LO 






So that 



F “ 




1 
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The coefficients of the sine terms in the Fourier series' for 
a square wave, magnitude .25, are tabulated in Column 4 of Table V. 
Each coefficient is identified with the harmonic to which it applies 
by the value of n in the same row, e.g. the coefficient of the 
sin 7 U-yt term, or seventh harmonic, is .0228. This procedure will 
be followed throughout the thesis, including the use of a second 
column for phase angle when required, and will be referred to as 
"tabulating the Fourier series." 

Having tabulated the Fourier series for the input square wave, 
and the frequency response of the linear system, we find the 
coefficients of the Fourier series for the output, °^(t), by 
multiplying each input coefficient by the corresponding* gain. 

The phase angles of the output series are, in this case, the phase 
shifts of the frequency response. The Fourier series for °^(t) 
is tabulated in Columns 5 an d 6 . From this tabulation we see, for 
instance, that the fifth harmonic term of the Fourier series for 
c< 1 (t) is .0775 Sin (5 u> t - 88.1°). 

The evaluation of (t), for 39 values of time up to 95 
percent of the half-period of the square wave, is tabulated in 
Table VI. The evaluation was made, as were all of those in this 
report, using the Wass, Hayman method. ^(t) is plotted in 
Fig. 6 together with the corresponding solution obtained by 
Brown-Edwards on a Short Analog Computer (see Ref. l). The 
oscillation is seen to be of high amplitude, due to the low damping. 
It is perhaps surprising that the approximation is as good as it is, 
because the VJass, Hayman method assumes that the oscillation will be 
nearly completely damped out by the end of the half cycle of the 
square wave, and it clearly is not. 

Solving Eq. 2-8 for [vj | = .25 we find 



w 



s 



.195 
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By Eq. 2-10, 



R " (.25) (.853) “ * 91? 

By Eq. 2-11, the evaluation of ^^(t) is found from the evaluation 
of °^(t). This is tabulated in Table VI, and plotted, with the 
analog computer solution of the non-linear equation, in Pig. 7. 

o< , is calculated by the method of Table IV. This calculation 
is shown in Table VII, and the resulting Fourier series has been 
transferred to Columns 7 and 8 of Table V. Multiplying the 
magnitudes by 

= (.917) 3 (.25) 3 = .01206 

and adding 180 ° to the phase angles gives the series for 
- \ y\ \ y or - °^2 9 this is tabulated in Columns 9 

and 10 . The frequency response ^/T(iuj), which in this case only 
differs from U(i^)/T(ito) by a constant multiplier, is tabulated 
in Columns 11 and 12 . By Eq. 2 - 16 , °< is found by applying the 

jO ^ ^ 

frequency response P/T(i«i) to the Fourier series for - , 

and the resulting Fourier series is tabulated in Columns 13 and 14. 
This is evaluated for 27 values of time in Table VI, and, by 
Eq. 2-17, combined with °^(t) to give °^^(t). °^(t) is 

plotted in Fig. 7. 

In Fig. 7 we see the relationship between the "reduced" linear 
solution ^^(t), the resu lt of the first iteration ^^(t), and the 
computer solution oc'(t). Superficially ^ doesn't look much 
better than particularly with respect to the amplitudes of 

the peaks and troughs. The times at which the overshoots occur 
have been brought into line however, and it will be shown in Section 3 
that the next iteration gives rapid convergence. The solution to 
the present example was not taken further because the digital computer 
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solution of a similar equation became available to the author, and 
the analog computer solution was thought not to have sufficient 
accuracy to justify its use as a standard with which to compare 
more refined approximations. 

SECTION 3: COMPARISON OF RESULTS WITH SOLUTION OBTAINED VJTTH 

A DIGITAL COMPUTER 



Mr. P. A. T. Christopher obtained solutions, from a Ferranti 
Mercury digital computer, of the equation: 

D, 2 x + a^Dx + a 0 x + bx^ = Q (3—1 ) 

where a^ = .6, 

a 0 = 1 ’ 

Q is a step function, 

Initial values are zero. 

These solutions were obtained for various values of b and magnitudes 
of Q. The accuracy of the solutions was such that, for present 
purposes, they may be considered exact. 

The output x(t), of Eq. 3-1, is dependent on four parameters: 
a l, a Q , b, and IqJ , the magnitude of the step input. It is shown 
in Appendix IV that by linear transformations of the variables x 
and t , the output can be made a function of only two parameters: 
the damping ratio % , and the non-dimensional input magnitude | Q q J . 
Thus the solution of Eq. 1-22, which is of the same form when 
7=0, could be found from the solution to Eq. 3-1, if V and 
| Q | for each of the equations were identical. Now the solutions 
to Eq. 3-1 include a range of values of | Q q | , but only one value 
of , i.e. .3, while the value of % for the missile at sea level 
is .175. Accepting this deficiency, we will apply the approximation 
method to Eq. 3-1 for the computer value of nearest that of 
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the case considered in Section 2. The results will not be 
transformed to the amplitude and time scales of the missile problem, 
because having already seen the approximate missile solution in 
Section 2, we are now only interested in the convergence of the 
approximations - which does not depend upon scale. 

The equivalent symbols in Eqs. 3-1 and 1-22 are tabulated 
below ( Z = o) : 



Eq. 3-1 



Eq. 1-22 






P> 

K 1 1 



Thus the non-dimensional input in terms relating to Eq. 1-22 is: 



$Jb (K x uji K ffi 1 

O 3/p + -7 

a o ( u>„2 ) 3/ 2 ^ 

The value of |Q q | corresponding to the case considered in Section 2, 
where the constants are those given in Table I for sea level flight 
and | | = .25, is: 

| Q ^ = l-8^J. l q346}.L..25l = >325 

° (24.4) 

The computer solutions nearest this are: 

(l) b 

Ul 

l Q ol 



= .5 

= .4 

= -4 fZ = -283 

d) 5/ 2 
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(2) b =1 

IQI = .4 

|Q \ = -4 fl = .400 

aW 2 

We choose the higher value, case (2), and proceed with the application 
of the approximation method to this case. In order to use the 
equations of Section 2 unaltered, we put Eq. 3-1 in the form of 
Eq. 2-2 by letting 

©4 = x 

= b = 1 

T(D) = D 2 + a,D + a = D 2 + .6D + 1 
U(D) = a Q = 1 




It I = IqI = .4 

The gain of the frequency response 

2 ^ , <„) 

T(iul) 1 - uj‘ + .6 uJi 

has a peak at CO = 1 , so we take oJL, - .2. Values of gain and 
phase shift of (3-2) are tabulated in Columns 2 and 3 of Table VIII, 
for the required values of nu3,. The Fourier series for a square 
wave, magnitude .4, is tabulated in Column 4. By Eq. 2-5 we 
find the Fourier series for (t), and this is tabulated in Columns 
5 and 6. The evaluation of for 22 values of time is 

tabulated in Table IX, and is plotted in Fig. 8 with the corresponding 
computer solution. The calculated solution appears to be a 

good approximation to the true output of txie linear system. 
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The steady-state output is found from Eq. 2-2 by setting the 
derivative terms equal to zero : 

(<*) 3 + (**>! - i-ii 

For hi- .4 this equation yields 

(<>O s = .355 

By Eq. 2-10, 

R - yifuT - • 838 

The evaluation of (t) is found from the evaluation of ° ^ (t) 
according to Eq. 2-11. This is tabulated in Table IX, and is 
plotted in Fig. 9* 

y 3 

The calculation of is performed in Table X. The 

lu 

resulting Fourier series is entered in Columns 7 and 8 of Table VIII. 
Multiplying the magnitudes by 

R 3 l*7l 3 = (.888) 3 (.4 ) 3 = .0447 

and adding 180° to the phase angles gives the series for 

3 3 v 3 rz 

-iv | h I or - , and this is tabulated in Columns 9 

K i u 2 

and 10. The frequency response O t is the same as 

U(i u))/T(iu) ) in this case because ft = u(it4->J = 1. This has 
already been tabulated in Columns 2 and 3* Thus is found 

according to Eq. 2-16, tabulated in Columns 11 and 12, and 
evaluated in Table IX. is obtained using Eq. 2-17, and plotted 

in Fig. 9. 

The next approximation requires a numerical harmonic analysis. 

3 3 

In Table IX, °( 0 and are cubed, and the difference ( °^ 0 - o< i) 

15 3 5 ^ 3 

is found for each of the 22 values of time. (o< ^ ) is 

plotted in Fig. 10. 
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Now harmonic analysis requires that the function being analysed 
be known throughout one period of the fundamental frequency uj^. 

So far only the first hal" 



Table XI. In this tabulation, the time t ^ is the time measured 



with those in the first half-cycle, and set them all equal to zero 
in the numerical harmonic analysis. This result is not surprising, 
because the oscillation in the second half-cycle takes place about 
the steady-state value of zero. Therefore the effect of the 
non-linearity is small, and the oscillation of the non-linear system 
is nearly the same as that of the linear one. This can probably be 

7 7 

taken as a general result, so that (ck - cxf J ) need not be 
evaluated for the second half-cycle. 

A 24 ordinate harmonic analysis of ( performed 

using the layout recommended by Wylie in Ref. 3, is shown in 
Table XII. The positions at which the ordinates were taken is 
shown by the vertical lines in Pig. 10. They were found by 
dividing the first half-cycle into 12 equal time intervals. 

The ordinates are labeled y^, y , etc., in accordance with Wylie*s 
notation. The ordinates were all multiplied by 1000 to make the 
harmonic analysis more manageable. The result of the analysis is 
a series in both sine and cosine terms which includes the 12th 
harmonic. Retaining only up to the 9th harmonic, the series was 
put into the form containing only sine terms, and this is tabulated 
in Columns 13 and 14 of Table VIII. is obtained by Eq. 2-26, 

and this is tabulated in Columns 15 and 16. It is evaluated in 
Table IX and, in accordance with Eq. 2-27, added to to give 

o<^. This is plotted in Pig. 9 with the other solutions. We 



evaluations of > and 
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see that the convergence is virtually complete in the region of the 
first peak, and probably good enough for most practical purposes in 
the entire range of the computer solution. 

This completes the demonstration of the method developed in 
Section 2. It would be desirable to find out how many iterations 
are necessary for larger values of input magnitude, and for what 
values of input magnitude the first iteration is sufficient, but 
this could not be done by the present author in the time available. 

SECTION 4: EQUATIONS OF COMPONENTS OF THE MISSILE CONTROL SYSTEM 

A. Incidence Vane 

One of the components of the control system, which will be 
designed in Section 5, is an incidence measuring vane at the nose 
of the missile. The design of a conical type vane suitable for 
this use is shown in Appendix V. The vane will be pivoted at its 
c.g. position. 

The relationship between the angle of attack of the missile 
and the output of the vane must be determined. Pig. 11 shows the 
angles and sign conventions used in the analysis. The angle ^ is 
measured by a transducer in the vane and is the output of the vane. 

It will be seen from the figure that no distinction is made 
between the velocity of the c.g. of the missile and the velocity of 
the vane. This is not strictly correct, of course, because the 
pitching velocity of the missile q creates a component of vane 
velocity perpendicular to the missile axis with magnitude 1^ q. 

It is easily shown, however, that for the configuration and flight 
conditions under consideration, this component of vane velocity has 
negligible effect on the vane velocity. For the example being 
considered, under steady-state conditions at an angle of attack of 
40°, the component of vane velocity due to pitching of the missile 
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is 5.3 ft/sec., and the difference between vane incidence and missile 
incidence is •12°. 

The equation of angular motion of the vane about its c.g. is: 



«V * \ ®v (4-1) 

The aerodynamic moment on the vane is taken as 

>V = C 1 + °2 * (4 - 2) 



where 

Cl = 

3°<v 

and & is the internal damping moment arising from the 
relative angular motion between the vane and the missile. It is 
assumed that a device within the vane exists to provide this damping. 

As in the case of the missile (see remarks following Eq, 1-3) 
the aerodynamic damping moment, that due to ^ is omitted on 
the assumption that it is small, but for the real reason that a 
value for it is not available. 

Combining Eqs. 4-1 and 4-2 gives: 

9 

Ci©< v + c 2 * = Iy (4-3) 



tf/e wish to reduce Eq. 4-3 to a linear equation in the two 
variables o< and % . From the figure, 



and 



= <*- i 



0 v = 



(4-4) 



y 



(4-5) 
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Substituting (4-4) and (4-5) in (4-3) gives 



C X {<* - t) + C 2 i = Iy (© M - K ) 



(4-6) 



Now © or q represents the pitching accel/eration of the missile, 

0 0 

and )( that of the vane relative to the missile. It is reasonable 



to assume, and it will be evident later, that the frequency of 
oscillation of the vane is much higher than that of the missile. 

This being the case, and providing that the magnitude of the missile 
oscillation is not excessive, then © ^ may be assumed negligible 
in comparison with V . Accordingly, Eq. 4-6 becomes 



( °< - )0 + c 2 v = - i t 



or 





(4-7) 



The constants and are now defined as follows: 

9 ^ 




2 



u) 



S' 2 = V 2 h “S 



Eq. 4-7 then becomes 





(4 -Q) 



The calculation of the constants will be found in Appendix VI. 
They are tabulated in Table I. 
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B. Actuator 

The equation for a typical elevator actuator and the constants 
of that equation were obtained from D. Morton, a student at the 
College of Aeronautics who had investigated the subject, 
equation is: 

•7 + 2 ^3 *\ = -C ^ 

b 3 3 

where ^ is the demanded elevator angle, 
is the elevator angle. 

The constants uJ^ and $ v are given i n Table I. 

SECTION 5 : DESIGN OF THE CLOSED-LOOP CONTROL SYSTEM 

A. The Linear System in Transfer Function Form 

The control system arrangement to be used is shown in Pig. 12 
in block diagram form. In this diagram: 

is the demanded missile angle of attack, the control 
system input, 

c< R is the reference input, equal to times 'K’j where 

the constant is the gain of the linear amplifier, 
o< is the missile incidence, 

is the demanded elevator angle, 

H is the elevator angle, 

]( is the output of the incidence measuring vane, 

and 

E = = Kj °< D - t (5-1) 

Equation 4-8 relates «< and V , Eq. 4-9 relates and , 



The 



D 



(4-9) 



and Eq. 1-22 relates o< and ^ . Dropping the non-linearity for 
the present so that the design problem is that of a linear system, 
Eq. 1-22 becomes 
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+ 2 5 1 



uJ, 



°< + ^ 



o< = 






( 1 + Z ) 



(5-2) 



The relation between E and ^ is, as yet, unspecified. 

The Laplace transformation, defined by: 

oO 

f(s) = J f(t) e~ St dt 

O 

where f(s) is the Laplace transform of f(t), will now be applied 
to Eqs. 4-8, 4-9, 5-1, and 5-2, assuming that initial values of 
all the variables and their derivatives are zero. Taking the 
Laplace transforms of these expressions - when initial values are 

zero - simply amounts to substituting s for d/dt and s for 
2 2 

d /dt . These equations then become, after rearranging, 



1(A 



Uj 



v(^) = 



ox 



* + 2 u)^ ^ + uj. 



'TU 



(5-3) 



r \U) 



=■ W {*') = 



2 



+ 2 5 



^ + tO 






(5-4) 



E(^) 



°< D ^) - 



(5-5) 



o< (/4-) _ 

Pf (U) 



K lP (^) = (l +^ 



Xt 2 + 2 £,**W -4- + 



(5-6) 
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where V(s), W(s), and P(s) are defined as indicated. ihe 

transform of the equation relating E and ^ ^ is given the 
form: 



In servomechanism terminology, V(s) is the transfer function of 
the vane. Similarly w(sj, r( s ), and C(s) are the 
transfer functions of the actuator, the missile aerodynamics, and 
the compensating network-amplifier. 

Following standard servomechanism technique, the control system 
may be represented by the block diagram in Fig. 13. 

From Eqs. 5-3, 5-4, 5-6, and 5-7 (or directly from Fig. 13): 



This is known as the open-loop transfer function, and the system 
comprised of the four elements whose transfer functions appear in it 
and which has an input E(t) and an output X(t) is the 
"open-loop system". 

Combining Eqs. 5-3, 5-5, and 5-8 to eliminate E(s) and Y(s) 
gives the closed-loop transfer function: 



= K, C(s) 

E(s) 2 



(5-7) 



where 



c(o) = 1 



Y 

3 (s) = K ± K 2 C(s) W(s) P(s) V(s) 



(5-8) 




(5-9) 



C(s) has not been specified. It has been included in order 
to provide the system with satisfactory stability characteristics. 
The stability of the system and the required form for C(s) will 
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be determined using frequency response technique. The choice of 
K ^ and will first be discussed. 

B. Steady-State Accuracy 



Suppose that ^ 



D 



is set at some steady value ( ) . The 

D o 



rules of Laplace transforms tell us that (in this type of system) 
takes on a steady value (<=><; after the transient motion 
has disappeared, and the ratio (°Og/(°<p)g is given by setting 
s = 0 in <=* Accordingly, from Eq. 5-9: 



ink. 



f-(0) 



K, K. 



2«3 



1 + K x K 2 



( 5 - 10 ) 



since 



c(o) = w(o) = P(o) = v(o) = 1 



Ideally the ratio of steady-state values should be one. 

That is the steady-state output should equal the steady-state 
demand. We see from Eq. 5-10 that this can be accomplished by 
a suitable choice of K^, regardless of the value of • 

But there is another important consideration: if the missile were 

to operate entirely at one flight condition, that is at a particular 
weight, c.g. position, speed, and altitude, then would always 
have the same value. But since missiles rarely operate under such 
conditions, we may assume that will not maintain a single 

value throughout the flight, and that will therefore 

vary (unless is made to automatically compensate it). It is 

evident in Eq. 5-10 that the higher the value of K ^ can be 

made, the less will be the variation in (,<=<) /(°< ) Q for a given 
change in K^. is determined by the aerodynamics of the 

missile, but may be chosen. 

Aside from considerations of amplifier design, there are two 
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factors which impose an upper limit on K^. One is that as 
is increased, the system becomes less stable, and finally 
unstable. So that the choice of must be a compromise between 

stability and steady-state accuracy. The other limit imposed upon 
is due to the presence of spurious inputs, or noise, which 
must not be amplified excessively. 

When has been chosen, may be chosen to give 



K 1 K 2*S = 1 

+ Kg 

for a given flight condition. 

If the steady-state equation for the aerodynamics, including 
the non-linear term, is combined with the steady-state equations 
for the other system components to eliminate all the variables 
except (°Og and > the following equation is obtained: 



(<=*) 

___S_ 

(c< ) " = 



(° < )s + 



ft 

(1 + k 1 k 2 ) 



(o<) s - 



K x K 2 
1 + K x K 2 



(o< ) 

K D J S 



( 5 - 11 ) 



We see that the larger is made, the smaller is the 

coefficient of the cubic term. Thus a high open loop gain, seen 
to be desirable in the linear system, reduces the effective 
non-linearity in the non-linear system. 



C. The Frequency Response of the Linear System 

The frequency response of a linear system at the frequency 
^0 is defined as the complex number obtained by substituting iu> for 
s in the transfer function relating the output to the input. It 
is usually expressed in polar form, with the magnitude in units of 
decibels. A number is expressed in decibels by multiplying its 
common logarithm by 20. 
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Now in a linear system a steady sinusoidal input of constant 
frequency cO will give rise to a steady sinusoidal output with 
the same frequency, after the transient motion has disappeared. 

The rules of Laplace transforms tell us that the ratio of the 
output sinusoidal wave amplitude to the input sinusoidal wave 
amplitude, or gain, is equal to the magnitude of the frequency 
response of the system at frequency uj ; and that the angle by 
which the output sinusoidal wave leads the input sinusoidal wave, 
or phase shift, is equal to the argument of the frequency response 
at frequency oo . Furthermore, by virtue of the properties of 
complex numbers, it can be shown that if the transfer function of 
a system is the product of two or more individual transfer functions, 
then the phase shift of the system is equal to the sum of the phase 
shifts of the individual transfer functions; and, if the gains are 
expressed in decibels, then the gain of the system is equal to the 
sum of the gains of the individual transfer functions. 

The frequency response curves of a system, sometimes called 
simply the frequency response, are graphs of the system gain and 
phase shift versus uj (or u> times some constant). 

It will be seen that the frequency response of a system may be 
found by substituting i cO for D in the differential equation of 
the system, since the transfer function is found by substituting 
s for D and the frequency response by substituting ioj for s. 
This procedure is used elsewhere in this thesis. 

As is shown in textbooks on servomechanisms (see Bibliography), 
the frequency response curves of the open-loop system determine the 
stability of the closed-loop system, and they may be used as the 
basis of trial-and-error design. 

Evaluating the individual frequency response curves and 
combining them in accordance with the rules previously stated, 
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V(i <o) W(iu)) was found and is plotted in Fig. 14. The 
parameter is used as the abscissa, and is plotted on a 

logarithmic scale for convenience. The constants are those 
given in Table I for flight at 50,000 f . It can be seen from 
Eq. 5-8 that P(i uo) V(iu>) W(i^) need only be multiplied by 
C(ito) to give the open-loop frequency response. It is 
now necessary to choose a particular compensating network. 

In ordinary servo design the choice of the compensator is 
based on the two main requirements: l) System stability, and 

2) Independence of the closed-loop system from changes in the 
constants of the components, i.e. high open-loop gain. But in 
the present example an additional requirement must be met: the 

system must have a frequency spread amenable to the Wass, Hayman 
method of transient response calculation (Appendix II ). In fact, 
it was at this point in the analysis that the writer first became 
aware of the frequency spread requirement, by having chosen a 
compensator that did not satisfy it. To demonstrate: This was 

a lead-lag network with the transfer function: 



C(s) = 



(1 



OJ 



*)(l 



^ I 



1.25 

'TXt 



O 



(l + 0(l + -125 *s ) 

Using the logmodang plot, Nichols chart technique (see Bibliography), 
the open-loop gain K ^ was set at 17 db, giving a closed-loop 

peak gain in the upper range of acceptable values - i.e. a lightly 
damped transient. was set at 1.14 db to make the closed-loop 

steady-state gain one. The closed-loop frequency response curves 
are plotted in Fig. 15. The V/ass, Hayman method was applied to 
this linear system to find the transient response for a unit step 
input of • The fundamental frequency of the input square wave 

was taken as .3 <*4*^ , one fifth of the frequency at the resonance 
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peak. The 13th harmonic was included. The transient response 
thus derived is plotted in Fig. 16. It appears to settle down to 
a steady-state value of approximately . 8 . But the actual 
steady-state value is one. The method has failed, because the 
frequency spread is too large. We see from Fig. 15 that the 
frequency spread, as defined in Appendix II, is approximately: 



6 



.08 



= 75 



Many other compensating networks were tried in an effort to 
obtain a satisfactory compromise between the three stated 
requirements. In fact, none was found, and it appears to the 
writer that the characteristics of the uncompensated open-loop 
system under consideration make the existence of such a compensator 
impossible. 



In order to continue the investigation of the non-linear 
problem, the actuator, vane and aerodynamic transfer functions were 
arbitrarily changed to make a satisfactory compromise possible. 

The actuator and vane transfer functions were set equal to one, and 
the aerodynamic damping ratio 5^ was changed from .0685 to .3. 
The block diagram of the simplified system is shown in Fig. 17. 

The compensator decided upon was a lead network-amplifier with 
transfer function; 

_^5 

1 + 5 



C(s) 



1 + 






cO, 



'*1 



s 



was set at 10 db, again giving a fairly high value of 
closed-loop peak gain, and requiring that = 2.38 db for unity 
closed-loop steady-state gain. The closed-loop frequency response 
is plotted in Fig. 18. It has a satisfactory frequency spread. 
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SECT 101m 6: THE TKANSIEHT RESPONSE OF THE MOK-LESak CLOSED- 

LOOP SYSTEM 



In the simplified system, that is for 
W(s) = V(s) = 1 



Eqs. 5-1 and 5-7 become, respectively, 

E(t) = TLj °< D (t) - «X(t) 



( 6 - 1 ) 



and 



(t) = h 2 C(D) E(t) (6-2) 

Equations 6-1 and 6-2 are applicable to the non-linear system as 
well as to the linear one. 

Using the same notation as in Section 2, the differential 
equation for the aerodynamics is: 

T(D) <=< = U(D) *| - (i J (6-5) 

Where 

T(D) = D 2 + 2 5 D + UCL (6-4) 

1 1 1 

U(D) = K x u>£ (1 + ZD) (6-5) 

Substituting (6-l) in (6-2) gives 

H (t) = K 2 C(D) [ K 3 o< D (t) - o<(t)] (6-6) 

Substituting (6-6) in (6-3), 



T(D) cK = K 2 C(D) U(D) [iCj °< D - o<^ 



(6-7) 
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This may be rearranged in the form: 

[t(d) + K 2 C(d) U(d)] o< = K 2 Kj C(d) U(d) °< d - ft (6-8) 

Now an approximation to the solution of Eq. 6-8 is the solution 
<=*^(t) of: 

[t(D) + K 2 C(D) U(D)] 2 = K2 C(D) U(d)* d ( 6 - 9 ) 



In transfer function form, Eq. 6-9 is 

^ Kp K_ C(s) U(s) 

_ 1 ( 5)= 2-2 

D T(s) + K 2 C(s) U(s) 

K 2 S C < s > iff) 

1 + *2 ^ wr 

Comparing Eqs. 6-4 and 6-5 with Eq. 5-6 we see that 



( 6 - 10 ) 



$t} - % 



Therefore : 



% K 2 C(s) P(s) 
1 + K 2 C(s) P(s) 



( 6 - 11 ) 



( 6 - 12 ) 



The right hand side of (6-12) is simply the closed-loop transfer 
function of the linear system (cf. Eq. 5-9) » so that °<^(t), the 
first approximation to the transient response of the non-linear 
system, is the response of the linear system - not an unexpected 
result. 
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Now by analogy with the argument presented in Section 2, we take 
the second approximation to be : 

o< z ( t) = R °q(t) (6-13) 

where 

R = V ; SN (6-14) 

^SL 

( c ^ / - is the steady-state output of the non- 

linear system, 

(0<) SL is the steady-state output of the linear 

system (equal to the magnitude of the input 
in the present example). 

The third approximation to o<vt) is °< .(t), the solution 
of: 

[t(d) + k 2 c(d) U(D)]«* 3 = K 2 c(d) u(d) °< d - ft °< 2 

Or, if °^..Ut) is the solution of: 

2 * 

[t(d) + K 2 C(d) U(d)]<* 3 * = - fi °<2 (6-16) 

then by superposition of the solutions of (6-9) and (6-16), 

= °<i + °< * (6-17) 

The solution of Eq. 6-16 follows the same lines as that of 
Eq. 2-16: we have o<^(t) in the form of a Fourier series: 

°h (t) = KJ \j- + tt 2 T sin (n ^ + ^ 

where K d I is the magnitude of the step input, and the A's 
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and B's represent the frequency response of the closed-loop system. 
As before we let 



©< M = / + 4 2 Sin < n V + V 

1 u 

so that 

%(t) = I»<bI ^ lo (t) 

Prom Eq. 6-13: 

°<' 2 (t) = R^t) = rI«< d | ^ (t) 

and 

°<2 = ^ I ^ dI 5 *1 



(6-18) 



We find the Fourier series for c*'. 



and then, by Eq. 6-18, the series for o( 
is found by applying the frequency response 



by the method of Section 2, 

3 



By Eq. 6-16, oC 



H(itd) = 



T(i«*>j + C(i<-o) U(i«->) 



3* 



(6-19) 



to the series for -<=*, . H(iuj) may be rearranged, using 

Eq. 6-11, to simplify its calculation. It becomes: 

K x K 2 Kj C(itO) P(i«*>) 

1 + c(iu>) P(ia>J \ 



H(i«0 ) = 



K, 



c(i»u ) U(i to. 



( 6 - 20 ) 



The bracketed term is the linear system closed-loop frequency 
response. 

The next approximation to o<(t), if it is desired to go 
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further, is ®^(t), the solution of: 

[ T(D) + K 2 c(D) UU>)] = K 2 Kj C(D) U(d)o< d - /3 o< (6-21) 

Adding and subtracting ft ‘K on the right hand side of Eq. b-21: 

[t(d) + K 2 c(d) u(d)] c< 4 = IC 2 k 3 c(d) U(i)) o<^ 

-P*? + p _ o0 (6-22) 

Let °^.( t ) be the solution of: 

[t(d) + K 2 C(D) U(d)] o< 4 „ = f 3 ( ^ - °< j ) (6-23) 

Then by superposition of the solutions of (6-15 ) and (6-23), 

c< 4 = + <* 4 * (6-24) 

3 3 

To solve Eq. 6-23, the Eourier series for ( °^ 2 - o<^ ) is 

found by numerical harmonic analysis in exactly the same way as 
outlined in Section 2 and demonstrated in Section 3* The 
frequency response H(ito) is applied to it to give °< 4# , and 
c< 4 is found by adding the evaluations of c<^ and o<^. 

Further iterations follow the same lines. 

Having calculated °^(t) and °<^(t) for a step function 
of any magnitude, it is easy to find their values for any other 
input magnitude. Let the subscript M a n denote the case which 
has been calculated and "b" the case which is to be calculated. 

The relationship between the linear solutions is, of course, 

o< (t) = ^b[ c< 1 (t) 

b l^al 3 



(6-25) 
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Prom Eq. 6-13 » 



c < 2 (t) 

a 



= E a °<1 M 



and. 



°<o 00 = 



R b “<,(*) 

b 



Substituting (6-25) in (6-27): 

<* 2 ( t ) = *1 M 

b l*®. I a 

Eliminating c<^ (t) between (6-26) and (6-28): 



<=<P (t) 

b 



Thus 



R bl 


I'XI 


R 

a 


H D al 


/M 




M 


[°<D 1 



o<2 0) 

a 






(6-26) 



(6-27) 



(6-28) 



(6-29) 



(6-30) 



Due to the linearity of Eq. 6-16, 

\ 

‘3\ 







(6-31) 



Having found o< (t) and C3 C.. (t), 
x b ^ b 

Eq. 6-17* No new Fourier evaluations have been performed. 



(t) follows from 

3b 



The only simplification possible in the next iteration is 



X 

that the evaluation of follows from Eq. 6-15 



3b 



must be found by cubing point by point, and a new harmonic 

analysis must be performed. 



We now apply this procedure to the simplified control 

system, first for a step input of magnitude .3 radian. 

The fundamental frequency ul is taken as .38 tcL one fifth 

* 1 

of the closed-loop frequency response peak in Fig. 18. The 
tabulations of Fourier series and frequency responses for this 
example will be found in Table XIII, and the evaluations in 
Table IXV. The Fourier series for a square wave, magnitude .3, 
is first tabulated, then the closed-loop frequency response. 

The Fourier series for the output of the linear system, (t), 
is found by applying the closed-loop frequency response to the 
square wave. We find the steady-state output of the non-linear 
system from Eq. 5-11, which, upon substitution of the constants, 
becomes : 



( °^SN + (°^SN “ 



( 6 - 



For 



find 



= K|= -3 



<*>SN = - 287 



Thus, from Eq. 6-14, 



.287 

.3 



.957 



°< 2 (t), found by Eq. 6-13, is evaluated. 

The Fourier series for , found by the method developed 

1 v 

in Section 2, is tabulated. 



R 5 |®< D | 3 = (.957) 3 {. 3 ? = .02365 

By Eq. 6-18, - is found from and tabulated. The 

frequency response H(iOL>) is tabulated, and the Fourier series 
for is obtained by Eq. 6-16, tabulated, and evaluated. 

By Eq. 6-17, the evaluation for (t) is found by adding those 
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for and 

1 3 

and °C(t) are plotted in Fig. 19. We note that 
<X^(t) is not very different from °< ,(t) and is therefore 
probably a close approximation to o<(t). The effect of the 
non-linearity is small, causing- a very slight increase in frequency. 

To illustrate a case where the effect of the non-linearity 
is large, we take = 1 radian. Although this is an 

unrealistic value of angle-of-attack, the result will be equivalent 
to that of a missile at low angle-of-attack but with a larger 
non-linearity. It is simply convenient to calculate it in this 
manner since the equations have been formulated on this basis. 

The evaluations are given in Table IXV. From Eq. 6-32 find 



Using the subscript "b" to denote this case, we have, from 
Eq. 6-14, 



The subscript "a" denotes the .3 input case, already evaluated. 






.757 



.757 



1 





(2.63) 5 = 18.19 



o^(t) is found by Eq. 6-31, and by Eq. 6-17. °^2^ an< * 

°<3b are plotted in Fig. 20. 
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One additional graph is presented in Pig. 21. Here the third 
approximations corresponding to both the .3 and the 1.0 input have 
been plotted, as has that for the .65 case, the calculations for 
which are not shown. These solutions have been ’’normalised" by 
dividing each ordinate by the corresponding steady-state value 
(°<)<^. 9 of course, is the same for all 

inputs, has also been plotted. This is a useful method for 
visualising the effect of increasing input magnitude on c^(t). 

It would also be useful for plotting other approximations, or the 
exact solutions. 

Improved approximations will not be demonstrated here, and 
the determination of the number of approximations necessary for 
acceptable convergence remains for the worker who thinks this method 
might be useful. 

Finally, it is interesting to note that a linear closed-loop 

system with two inputs can be described which is equivalent to the 

non-linear system. The block diagram of this system is shown in 

Fig. 22. Comparison with Eq. 6-8 will prove its equivalence. 

Fig. 22 is, in fact, the general outline of a possible analog 

3 

computer set-up where o< is found by putting o< through a cubing 
device. 
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DISCUSSION 



It has been mentioned that a quick method for determining the 
effect of non-linear stiffness was desired* Specifically, it was 
hoped that few terms would be required in the first iteration, and 
that one iteration would be sufficient for acceptable accuracy, at 
least in the first overshoot region. In fact, the number of 
significant terms was not small, and a single iteration was not 
sufficient no matter how many terms were retained, except for small 
inputs. 

The frequency spread requirement, discussed in Appendix II, 
places a limitation on the range of systems to which the Wass, 

Hayman method may be applied. And the Wass, Hayman method must be 
applicable for the present method to work. If the frequency spread 
is of the order of 20, the number of terms in the Fourier series 
could be increased to handle it, but if the original control system 
designed here is at all typical, the frequency spread is likely to 
be at least 50, requiring at least 25 terms in the Fourier series 
for the linear response, and more for the non-linear response. 

The need would have to be great to justify the effort of evaluating 
such series. 

In Section 5 an unrealistic value of aerody nami c damping was 
arbitrarily chosen to produce a control system to wnich the present 
method could be applied. It is possible that a more justifiable 
change in the control system could have been made, but this was not 
done due to a necessity to find an acceptable system quickly. It 
remains to be seen whether the present method can be used in any 
practical cases. 

Although only "hard" stiffness has been considered, i.e. 
positive /3 , the method is equally applicable to soft stiffness, 
or negative ft . In the latter case the presence of static 
instability at high angle-of-attack must be considered. 
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CONCLUSIONS 



1. Use of the iterative technique developed is 
capable of giving the transient response, to an acceptable degree 
of accuracy, if the frequency response of the system meets the 
frequency spread requirement. 

2. In an example, acceptable accuracy was obtained 
after two iterations. The number of iterations necessary when 
the non-linear effect is large was not determined. 

3. The labor required to evaluate the Fourier series 
and to perform the numerical harmonic analyses necessary in the 
second and subsequent iterations is considerable, thus limiting 
the practical value of the method. 
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APPENDIX I 



CALCULATION OF THE AERODYNAMIC CONSTANTS 

Brown-Edwards, in Ref. 1, found, for the missile under 
consideration and for a c.g. position 8 body diameters aft of the 
body apex, the following values: 



C L = 2.87 - .222 



( 1 - 1 ) 



C M - - 1.14 - 2.65 + .97 n 

c.g. 1 



Fixing the body diameter 1 at one foot, 
typical missile: 



rm. 



12.05 



. v 2 

lb. sec 
ft. 



he estimated for a 



( 1 - 2 ) 



and 

I = 250 lb. ft. sec. 2 

2 

The reference area S iS 9 ft . 

We take the case of flight at 50,000' at Each 3. The air 
density p is .000362 lb.sec 2 /ft^ and the speed of sound c 
is 968 ft/sec. Thus the velocity is: 



V = 3 (968) 



2900 ft/sec. 



The dynamic pressure is 



Q = TfV 2 



£(. 000362 ) ( 2900) 2 




and 



QS = (1520)(9) = 13,700 lb. 



Now 



7 _ 

~ 3« 



3C t 

= " $ s ^ = " (13 >700; (2. 87) = - 39,400 lb. 



where '9c./3o< is obtained from Eq. 1-1. 
Jj 



9L 

d \ 



dC 

2 ~ ^ = -03--= - U3,700)(- .222) = ;>,040 lb. 

d>i o»| 



where 'dCj/'drj is obtained from Eq. 1-1 

M = pi = qs£ = (13,700)(l)(.97) = 13,300 lb. ft. 

v o\ o ^ 

where 3 C M /0 ^ is obtained from Eq. 1-2. The moment due 

c.g. ' 

to o< is, from Eq. 1-2, 

M"(©<) = Q <5 £ (-1.14 - 2.65 c< 5 ) 



= (13,700)(l)(- 1.14<* - 2.65oO 



= (-15,600 °< - 36,3OO<*0 lb. ft. 



Comparing this with Eq. 1-14, we see that 
M* = -15,600 lb. ft. 

and 

M x = -36,300 lb. ft. 

Using Eqs. 1-17 to 1-21, 






*1? - 15,600 = 67.8 sec 2 . 

T“ 250 



39,400 

2 (12.05) (2900) (8. 23) 



U> '^ 1 = 8.23 sec 



-1 



5i * - 



2 -"m-V 



= .0685 



K, = 



i-i 



lu) 2 

•x 

1 



13,300 



(230)(67.8) 



= .853 



z = 



3,040 






(12. 05) (2900) (.853) (67.8) 



ft = 



.001505 sec. 
M o<3 



I 



36,300 

230 



158 sec 



-2 



APPENDIX II 



THE WASS. HaYIlAH METHOD OP DERIVING THE TRANSIENT RESPONSE OF 
A LINEAR SYSTEM FROM THE FREQUENCY RESPONSE 

Reference 2 presents a method of deriving the approximate 
response of a linear system to a step input from the frequency 
response. The method is to find the transient response to a 
square wave input instead of to a step input. The justification 

for such a manouver is best explained with reference to Fig. 23 . 
Figure 23a is the graph of a step input, and Fig. 23b is the 
corresponding output of a typical system. Figure 23c is a square 
wave input of period T, and Fig. 23d is the corresponding output. 

It is evident that if the half -period T/2 is greater than the 
time required for the transient to completely damp out, then any 
of the half cycles of Fig. 23d is the same as the transient part 
of Fig. 23b. Of course in linear systems the transient never 
completely damps out. Therefore Fig. 23d is really an 
approximation to Fig. 23b, its accuracy depending on the extent 
to which the transient subsides in the time l'/2. 

The reason for using a square wave is that it can be expressed 
as a Fourier series. The Fourier series for a square wave F^(t) 
with equal mark-to-space ratio, minimum value zero, maximum value 
unity, and fundamental frequency tO , can be written: 

F.(t) = t+ £ 'S’ i sin n<t 

1 TT k__i n 

n = 1, 3, 5, . • . , 00 

(Wass and Hayman use the symbol n differently - this will be 
explained shortly). If Ih(t) is the input to a linear system 
the magnitude and phase of each component will be changed 
according to the frequency response of the system, and the 



output will be 



p (t) 

o J 



A o 2 X A n . , 

— + ^ Z — sm (n<o t + 



B ) 
rr 



n — 1 , p , 5 > ••• 



where A is the gain of the system at frequency n uj , 
n F 

B is the phase shift of the system at frequency nu> . 
n j? 

In order to use this method numerically, the summation must 
be limited to a reasonably low number of terms. Wass and Hayman 
state that it has been found by experience that a satisfactory 
compromise between labor and accuracy is achieved by considering 
up to the 11th, or sometimes 13th, harmonic of the fundamental 
square wave frequency. Having decided on the number of terms 
to be summed, a value for the fundamental frequency must be 

si 

chosen. Wass and dayman suggest that for systems whose gain 
curves exhibit a definite peak, to , should be one fifth of the 
frequency of the peak. Suggestions are also made with respect 
to other types of systems. 

v*ass and Hayman do not point out that taking the summations 
up to even the 13th harmonic will fail to give satisfactory 
accuracy in certain types of systems. They do infer that the 
transient response is mainly determined by the characteristics of 
the frequency response at frequencies lower than that at which 
the gain is 15 db below the zero frequency gain. They do not 
point out that the gain and phase shift at the fundamental 
frequency must not differ greatly from their zero frequency values. 
For conciseness we may use the term "frequency spread" to denote 
the ratio of the -15 db frequency to the frequency at which the 
gain or phase first differs appreciably from its zero frequency 
value. Thus the Wass, Hayman method, summing up to the 13th 



harmonic, is not applicable to systems whose frequency spread is 
greater than 13* 

An important part of the Wass, Hayman report is a table, the 
M ^ 11 Table, which may be used to simplify the evaluation of Fourier 
series. This is a table giving the values of n<-o t* in degrees, 
for 39 values of non-dimensional time t*, and for odd values of 
n from 1 to 13. Even values of n are not included because the 
Fourier series for the responses of linear systems to square wave 
inputs contain only odd harmonics. Since the evaluation of even 
harmonics is required in this thesis, an extension to the if Table 
is required, and this is given, for even values of n from 2 to 8, 
in Table XV. Wass and Hayman* s use of the symbol n precludes 
the possibility of even harmonics, and this is the reason it is 
used differently here. 

bass and Hayman present a design for a circular computer which 
reduces considerably the time required to evaluate a Fourier 
series. The computer suggested consists of a circular base with 
various lines and scales marked on it, two arms which pivot at the 
center of the circular base, and a cursor which slides along one 
of the arms. The author constructed and used such a computer, and 
eventually devised a different design which was superior in several 
respects. Fig. 24 is a photograph of this computer. It was 
originally a "Wind and Navigation Calculator Mk. I", and all the 
scales around the periphery are superfluous except tue 360 degree 
one. It consists of an opaque base and a transparent top, 
pivoted at the center. The horizontal lines, the numbers running 
down the center, and the arrowhead are etched on the base. The 
360° scale is etched on the top. The top is roughened so that it 

will take pencil marks. To evaluate the term 

A = B sin (nuj^t * + © ) 



the procedure is as follows: 



1. Rotate the top until the arrowhead points to 

(90° - © ) on the 360° scale. 

2. Make a pencil mark, on the top, at the point 
where the value of B appears on one of the vertical scales 
above the center. 

3. Rotate the top until the arrowhead points to 
nco^t 1 on the 360° scale. 

4 . Follow the new position of the pencil mark along 
the horizontal lines to the same scale used to place the mark, 
and read A. It is positive if above center, negative if below. 

5. Repeat steps 3 and 4 for other values of n^t'. 

I 

The new computer is superior to the Wass, Hayman version in 
these respects: 

a. The two arms of the Wass, Hayman computer must be turned 
through several revolutions without moving relative to each other, 
yet the angle between them must be adjustable. A design 
incorporating such a feature may be difficult to achieve, particularly 
for the home builder. The new computer obviates this difficulty. 

b. The angle#* between the arms and the position of the 
cursor must be reset for each new term of the Fourier series. 

With the new computer a mark for each tercn of the series can be 
made at the outset, each one being labeled appropriately. This 
may not seem important, but in practice it is. Mistakes can be 
corrected easily because the setting is not "erased" for the next 
term as with the Wass, Hayman version, and evaluations can be made 
for additional values of time, without resetting, after several 
points have been plotted. 



APPENDIX III 



GRAPHICAL ADD Ii ION 09 nIXE riABI'IONICo 
If, in the expression 

e sin ( ^>t + ©) + f sin («->t + <p) + g sin {u>t + if* ) 

e, f, g, ©, <p, and y are known constants, then an equivalent 

expression is 

h sin (u) t + ) 

where h and ^ are determined according to the construction 
below. 




This method is valid for any number of terms in the original 
expression. 



APPENDIX IV 



TRANSFORMATION OF THE CUBIC EQUATION 



Given the equation: 



d u w 3 

dt 2 + a l c it + & o* + b * 



Make the substitutions: 



t = fa 



z 



and 









Then Eq. 1-1 becomes: 



* _!l jLs + +> 



We shall call the quantity 

Q fb 



Q. 



o - 



7^2 

o 



fa 

1 n 



= Q (I- 



Qfb 

~ a 2/ 2 
0 



the non-dimensional input, 
twice the damping ratio 5 



The coefficient 



is equal to 



APPENDIX V 



DESIGN OF THE INCIDENCE VANE 



The vane is conical, and it is necessary that the c.g. be as 
far forward as possible. The configuration to be used is that 
shown in Fig. 25. The nose will be solid steel, and the base will 
be a thin shell of light-weight plastic. The optimum ratio of 
nose length to tail length, for most forward c.g. position, is that 
which places the c.g. at the junction of the nose and tail. The 
equation for this condition was found to be: 



_t 

L 



.1 sin O 

4 /°2 (2 n 4 - 3 n 3 + n) 



where 



t 


is 


the 


L 


is 


the 


A 


is 


the 


(* 2 


is 


the 


© 


is 


the 


L 

n 


is 


the 



base thickness, 
cone length, 

density of the nose material, 
density of the base material, 
cone semi -angle, 

distance of the c.g. from the 



apex. 



The density of steel is about .3 lb/cu in., and that of a high- 
temperature silicone plastic is about .05 lb/cu in. (p. 66 of 
Ref. 4)* Using these values, the compromise chosen was: 



© = 22.5° 
n = 3 
L = 3" 

t = .0614" 

It is questionable whether the base is thick enough for adequate 
strength, but it will be assumed that such a design is practicable. 



APPENDIX VI 



CALCULATION Or* TriE VaKS CffiSTA LTS 

Ref, 5 gives experimental pitching moment data for cones in 
supersonic flow in the form of graphs of versus , where 
the pitching moment about the nose is 

^OSE = U Q r b 1 c C L 

and 

r^ = base radius, 
l c = cone length. 

The curves of versus were virtually linear for the conditions 
applicable to tne present case, and the value of '9C T /3°< about 
the nose was .020. The center of pressure remained very close to a 
position 2/3 1 back from the nose, so that the moment curve slope 

for a cone hinged at 1/3 1 would be half of .020, or .010. 

From Appendix I: 

w = 1520 lb/ft 2 

Erom Aj.wena.ix v 

l c = .25 ft. 

A f = .0107 ft 2 , 

b 




= - TTU520j(.0107)(.25)(.0l) 

= - .128 ft. lb. 

p = 1.70 x 10 ^ ft. lb. sec? 

(Eound by integration, calculation not shown). 



.128 



iO 



n„ 



_i = 



= 7,530 sec 



-2 



1.70 



10 



-5 



= 86.9 sec 

n 2 

5 2 = .5 

(Chosen arbitrarily on the assumption that the vane damping device 
can be so designed.) 

C 2 = 2 IyU> = 2(1.70 x 10~ 5 )(86.9)(. 5) 

= .00146 ft. lb. sec. 



t 
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